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Quantum fluctuations in time-dependent, harmonically-trapped Bose-Einstein condensates are 
studied within Bogoliubov theory. An eigenmode expansion of the linear field operators permits 
the diagonalization of the Bogoliubov-de Gennes equation for a stationary condensate. When trap 
frequency or interaction strength are varied, the inhomogeneity of the background gives rise to off- 
diagonal coupling terms between different modes. This coupling is negligible for low energies, i.e., 
in the hydrodynamic regime, and an effective space-time metric can be introduced. The influence 
of the inter-mode coupling will be demonstrated in an example, where I calculate the quasi-particle 
number for a quasi-one-dimensional Bose-Einstein condensate subject to an exponential sweep of 
interaction strength and trap frequency. 



I. INTRODUCTION 

Ultracold atomic gases offer various opportunities for 
the study of interacting many-body quantum systems 
in a well-controlled environment d 0|. For instance, 
the Bose-Hubbard model - a simplified description for 
bosons in a periodic potential - can be studied with 
Bose-Einstein condensates confined in optical lattices 
0,3- Quantum gases have also gained much attention 
lately regardin g th e emergence of an effectiv e sp ace-time 
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their low-energy phase fluctuations obey the same co- 
variant field equations as a scalar quantum field in a cer- 
tain curved space-time. Hence, the study of phonons 
in this laboratory system might shed some light on as- 
pects of cosmic quantum effects, e.g., Hawking radia- 
tion [H, E3, El EE H3] or the freezing and amplifica- 
tion of quantum fluctuations in expanding spacetimes 
S S EE EI E3, El HI HI- Although the fluctuations 
in Bose-Einstein condensates are usually small, it has re- 
cently become possible in experiments to go beyond the 
classical ord er p arameter and resolve signatures of the 
fluctuations [IS M, US EE E3 • 

Theoretically, the fluctuations in a Bose-Einstein con- 
densate are usually treated as small perturbations of the 
mean field. The solution of the coupled field equations 
is rather demanding and often requires further approxi- 
mations, especially for time-dependent condensates. The 
Hartree-Fock-Bogoliubov method, see, e.g., [11], permits 
in principle the self-consistent propagation of the mean 
field and the quantum correlations for arbitrary time- 
dependences of the trap potential or interaction strength. 
But the scaling of the numerics with system size often 
limits the actual calculations to a low number of dimen- 
sions, certain symmetries, or a short time interval. Ther- 
mal condensates might be studied using the projected 
Gross-Pitaevskii equation [2§| , where the low-energy part 
of the fluctuations is expanded into an arbitrary basis 
and, in view of the large thermal occupation, treated 
classically; higher excitations as well as the vacuum effect 
are omitted. On the other hand, studies in the context of 
expanding spacetimes, e.g., [S S IS IS EE EI E2] , indeed 
focus on the quantum fluctuations but often assume a ho- 



mogeneous background or start with the hydrodynamic 
action, which is only valid on scales longer than the heal- 
ing length. 



In this article, I will discuss the evolution of the 
quantum fluctuations in trapped time-dependent Bose- 
Einstein condensates. The linear field operators will be 
expanded into their eigenmodes thus permitting the diag- 
onalization of the (initial) evolution equations. Although 
basis expansions of the field operator are frequently used 
(e.g., in [IS these references usually consider the 

harmonic oscillator eigenfunctions, a large number of 
which must be used in order to describe the excitations 
properly. By adopting the eigenmodes, a much smaller 
part of the basis needs to be considered and many more 
situations will become numerically feasible. (Note, how- 
ever, that in order to obtain the fluctuation eigenmodes, 
a relatively large number of oscillator functions must be 
employed - but they need not be propagated.) 



This Article is organized as follows. Section [TT] re- 
views the field equations and their linearization. The 
Bogoliubov-de Gennes equation for the linear quantum 
fluctuations can be diagonalized by an eigenmode ex- 
pansion, which will be performed in Sec. Mil How- 
ever, as soon as trap frequency or interaction strength 
are varied, off-diagonal terms appear. This coupling of 
different modes is negligible for excitations with ener- 
gies much smaller than the chemical potential even in 
time-dependent condensates, as will be shown in Sec. 
\TV\ where the order parameter is treated in the Thomas- 
Fermi approximation. It is then also possible to establish 
the analogy between phase fluctuations and a massless 
scalar field in a certain curved spacetime. If the excita- 
tion energies are of the same order as the chemical po- 
tential, the coupling of different modes might lead to a 
population transfer, which will be illustrated in an exam- 
ple in Sec. El 
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II. FIELD EQUATIONS 
A. Scaling transformation 

In dimensionless units, the field operator $ of a 
trapped (quasi)-Z?-dimensional Bose-Einstein condensate 
obeys the non-linear Schrodinger field equation [1H 
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^-+c 2 (t)^+ 5 (t)*t^ 



(1) 



with D-dimensional coupling strength g(t) oc a s . By 
changing the s-wave scattering length a s through Fes- 
hbach resonances [32j or by varying the trap frequency 
Lo(t), an external time-dependence can be prescribed 
on the condensate. The gas cloud will adapt to these 
changes and it will either expand or contract and with it 
the quasi-particle excitations residing upon it. A part of 
this background motion can be accounted for by trans- 
forming to new coordinates x = r/b(t) with scale factor 
b(t). The field operator then reads [331 ] 



b D/2 



(2) 



The phase $ = (r 2 /2)b/b is chosen such as to gener- 
ate an isotropic velocity field V$ = rb/b which (at least 
partially) describes the expansion/contraction of the con- 
densate. If the scale factor b(t) obeys 
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dt 2 
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g(t) b 2-D 
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(3) 



with go being the initial value of the coupling strength, 
a scaled field equation follows 



+ f 



*l>, (4) 



where trapping and interaction terms have acquired the 
same time-dependent pre-factor f 2 (t) and all other coef- 
ficients are time-independent. (The scale factor b(t) on 
the left hand side might be included into a redefined time 
dr = dt/b 2 , see Sec. llVBl ) 



B. Linearization 

For large particle numbers N, one might formally ex- 
pand the field operator into inverse powers of y/N [34[ 
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+ x + C 
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N 



(5) 



Here, A and N = A' A are the atomic operators. They 
commute with the linear \ = O(N ) and higher-order 
quantum excitations C and thus yield the exact conser- 
vation of particle number. The order parameter ipo = 



0(yN) in the center of the trap but diminishes towards 
the edge of the condensate. Insertion of the expansion © 
into the scaled Heisenberg equation (@]) yields the Gross- 
Pitaevskii equation for the classical background tpo [35| 



, 2 a 



<,'o 



(6) 



The linear quantum fluctuations \ obey the Bogoliubov- 
de Gennes equation [36] 



-^f + / 2 (y + 2 5o |^oI 2 



X + / 2 5o'0 2 X t , 
(7) 



and the residual terms comprise the equation of motion 
for C 
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V 2 dt C 
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2/2 ' 2 

+ g (2^oX f X + r X 2 + X f X 2 ) + O(go() . 



(«) 



These higher orders C must remain small in order for 
the mean-field expansion J5|) to be valid, i.e., for the lin- 
earized equation |(7|) to be applicable. This means that 
the terms involving products of x must remain small be- 
cause they act as source terms for higher orders (. 

From the Gross-Pitaevskii equation ([6]), I can infer 
when the evolution of the order parameter ipo is solely 
described by the scale factor b(t): apart from the trivial 
case f 2 = const, this occurs only when the spatial deriva- 
tives can be neglected with respect to the interaction and 
trapping terms, V 2 ^ -C f 2 (x 2 + 2go\ipo\ 2 )ipo, i.e., in the 
Thomas-Fermi approximation. Density and phase of the 
order parameter tfjo — J~qo then assume the form 



TF 

Qo 



Ho ~ x 2 /2 



,9o 



e(w> - x 2 /2) 



dt 



, f 2 (f) 
b 2 {t>) 



(9) 



where the Heaviside step function 0(/io — x 2 /2) is 1 for 
> x 2 /2 and elsewhere. In this approximation, the 
motion of the classical background becomes stationary 
and the scaled coordinates x are co-moving with the con- 
densate. 

The Bogoliubov-de Gennes equation can be tackled 
by introducing self-adjoint operators 



X+ = e-^x + e l *°x\ 



X- 



2i 



(10) 



with 4>o = argV'o being the phase of the order param- 
eter. These operators resemble (relative) density and 
phase fluctuations Sg/go = X+/\fQo an d 5<t> = X-/\/Qo 
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up to the prefactor l/^/go. Since this prefactor even- 
tually becomes large near (and beyond) the surface of 
the condensate, the smallness of 5g/ qq and S(f> cannot be 
ensured. Therefore, I will stick to x± in the following, 
but still refer to them as density and phase fluctuations. 
They obey 



a 



1 



V— +v V x + ~{V x v 



d 



1 



X- = -K-+X+ , 
X+ =IC-X-, (11) 



where the velocity field vq = V x </>o results from the resid- 
ual background phase beyond the Thomas- Fermi approx- 
imation J9]). However, Vq is small for go > and if / 
does not change too swiftly, because an almost homo- 
geneous phase will develop with only small deviations 
near the boundary of the condensate. Whereas attrac- 
tive go < invalidate the Thomas-Fermi approximation 
and generally v ^= even in the center of the trap. The 
differential operators on the right-hand sides 



K+ = 
K.- = 



V X 

2 



x 

T 

_VJ , «o , .a 
2 2 1 V 2 
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(12) 



generally do not commute for inhomogeneous conden- 
sates V x Qo 7^ 

[£+,£_] = fgo {(V^ ) + 2(V xQo )V x } . (13) 



III. EIGENMODE EXPANSION 

In order to define the initial state unambiguously, I 
will assume that the condensate is at rest before t = t- m . 
Then b = 1, b = 0, f 2 = 1, and v = such that the left 
hand sides of Eqs. (fTT]l reduce to partial time derivatives 
and the initial eigenmode equations for x± follow 



dt 2 ' 
dt 2 



X- 



-K.-K.+X+ , 
-IC+IC^x- ■ 



(14) 



Because K.+JC- ^ JC-JC+ for inhomogeneous conden- 
sates, cf. Eq. lfl3|). density and phase fluctuations of each 
mode must have different space dependences. This leads 
to the expansions [13] (I will adopt the sum convention 
throughout this Article for brevity; any indices appearing 
only on one side of the equation are to be summed) 



X+(x,t) = h+(x)X+(t), 
X-(x,t) = h-(x)X-(t) 



(15) 



of X± into different eigenmode bases {h+} and {h n }, 
see Appendix [A] for more details on how to obtain . 



Usually, these two bases are neither orthogonal nor nor- 
malized, J h+h^ ^ 5 nm 7^ f h-n^mi but instead can be 
chosen to be dual to each other 



d D x h^(x)h m (x) = 5, n 



(16) 



Note that this condition does not fix the norm of but 
still permits the multiplication by an arbitrary factor, 
/i+ — > A n /i+ and h~ — > (1/ K n )h~. Observables must be 
unaffected by this ambiguity, see App. iBl 

Insertion of the eigenmode expansion (flB"! into the lin- 
ear field equations ifTTj) yields a set of coupled first-order 
differential equations 



d_ 

dt' 



b di X « 



2B inn 

(t)X-+Vmn(t)X+. 



(17) 



with time-dependent coefficients. The symmetric matri- 
ces 



Anm (^) 
$nm (t) 



d D x KllC+h 



d xh„ /C_ h„ 



(18) 



are initially diagonal A n m{tin) = A n (ti n )5 nm and 
Bnm(tm) = B n (ti n )5 nrn . At later times, they acquire off- 
diagonal elements because of the different-time commu- 
tators [K+{t),K + {t')] ^ and [JC_(t),JC_(t')] + when 
JC±{t) ^ JC±{t'). The velocity coupling matrix 



d D xht 



Vo^x 



(19) 



is not symmetric, but vanishes for homogeneous phases 
0o of the order parameter, e.g., initially or in the 
Thomas-Fermi approximation (J9j) . For repulsive go > 
and slow variations of interaction strength g(t) and trap 
frequency ui(t), the order parameter phase is homoge- 
neous except for small ripples near the boundary of the 
condensate such that the matrix V nm is usually negligi- 
ble. 

From the evolution equations lfl7| with the initially 
diagonal coupling matrices (fl8j) . the introduction of 
bosonic operators (v n and a n creating or annihilating an 
initial quasi-particle is straightforward 



X-(t) = F™(t)a m + F™(t)c 
X+(t) = GZ\t)a m +G™{t)a, 



(20) 



Here, a bar shall denote complex conjugation, e.g., F£ 
(F™)*. The coefficients obey the initial values 



A,. 



2a 



" $nm i C n (^in) ^ 



2A r 



(21) 



where the phase has been appropriately chosen and the 
frequencies Vl m — ^/A m B m . The upper index of the co- 
efficients F™ and G™ labels the mode, while the lower 
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index denotes the components of this particular mode 
when expanded in a certain basis, e.g., the initial eigen- 
functions {h„}. 

Since the coupling matrices lfl~8|) become non-diagonal 
even for slow (adiabatically) variation of the trap fre- 
quency oj(t) or coupling strength g(t), the initial bases 
{h„ } cannot represent the eigenmodes at later times. Al- 
though /i„ might be employed in order to calculate the 
spatial correlation functions, see Appendix the use of 
these functions might be misleading regarding the cor- 
relations between different modes. Furthermore, when 
probing the excitations using, e.g., the scheme proposed 
in [39], the proper particles defined at the time of mea- 
surement will be detected and not the initial ones. 

The particle definition in time-dependent background 
is a non-trivial task, see, e.g., [22]. Nonetheless, it is 
always possible to expand x± 



X±(x,t) = h% ti (x)X± tl (t) 



(22) 



into bases {h„. tl }, which are defined such that 
the coupling matrices Anm-M = / dPx h^. ti IC+h^ n . ti and 
Bnm-tx = J d D xh~. ti K.-h m . ti , cf. fl8|) . become diagonal 
at any particular instant t\ 

Anm\t\(tl) — A n -i 1 5 nm , ^nm;fi(^l) — B n -i x fa nTa ■ (23) 

Of course, the velocity term Vnm-.t! is then generally 
non-diagonal and the evolution equations of the differ- 
ent modes will not exactly decouple at this particular 
instant t\. But one should bear in mind that measure- 
ment occurs usually in an adiabatic region, where the 
external parameters are only slowly-varying functions of 
time. Then, the background phase is approximately ho- 
mogeneous and the velocity t>o ~ 0. Hence, quasi-particle 
creators and annihilators might be introduced analogous 
to Eq. (HI 



(24) 



where ft n;tl = y/A n . tl B n . tl . 



IV. THOMAS-FERMI APPROXIMATION AND 
EFFECTIVE SPACETIME 

In the previous sections, I made no approximations ex- 
cept for the linearization <[5j) and the assumption of an 
isotropic trap. The formalism is, in principle, applicable 
for arbitrary variations of trap frequency ui(t) and inter- 
actions g(t). To this end, it would be necessary to solve 
the Gross-Pitaevskii equation §§§ and the linear evolution 
equations lfl~7|) simultaneously. The numerical solution is 
complicated by the fact that the coupling matrices lfl~8|) 



and l|T9|) need to be calculated at each time step. Some of 
the numerical difficulties can be circumvent by adopting 
the Thomas-Fermi profile , where density and phase of 
the background become time-independent (in the coordi- 
nates x) and thus require the calculation of the coupling 
matrices only once. Despite some shortcomings regarding 
the dynamics of the order parameter, this approximation 
is usually applicable for repulsive interactions and in the 
center of the trap, but becomes inaccurate towards the 
surface of the condensate, where the quantum pressure 
°c (Vy/go) 2 is relevant. 



A. Coupled evolution equations 

Within the Thomas- Fermi approximation J9]), the co- 
ordinate transformation r — ► x associated with the scal- 
ing transformation <[2j) renders the background density 
time- independent gJ F = 0, while the phase becomes ho- 
mogeneous TF — and thus V^m = 0. The inte- 
grals of the coupling matrices (HU simplify considerably 
and the evolution equations can be cast into the form 

-26 2 ^A- = A n X+ + (f - l)M nm X+ , 

i& 2 ^A+ = B n X- + (f l)M nm X m , (25) 

i.e., Anm and B nm can be split into time- independent 
diagonal parts, cf. Eqs. (fl8|) . 



f ( V 2 x 2 \ 

A n 5 nm = / d D x h+ ( — + — + 3g Qo F - (j, J /i+ , 

f ( V 2 x 2 \ 

B n S„ m = / d D x h~i — y + Y + 9oQo F - Ho J Kn 

(26) 

and constant, symmetric, non-diagonal coupling matrices 

t „ . /cr 2 

M n in 
■A/rim 



J d D x h~ ( y + g gj ¥ - Ho) Kn (27) 

with time-dependent prefactors f 2 (t) — 1. The external 
variation of trap frequency u(t) and coupling strength 
g(t) is solely encoded in the scale factor bit) and the 
scalar function f 2 (t) = g(t)b 2 ~ D . Note also that the co- 
efficients |27j) and thus also the evolution equations lf25j) 
for the fluctuations are independent of the initial coupling 
strength go. The addend goQo F appearing in the paren- 
theses of Eqs. f27j) and (|26|) can be expressed through 
the chemical potential goQo F = {Ho ~ a; 2 /2)6(a;TF — x), 
where xtf = V%Ho, cf. the Thomas- Fermi equation ([9]). 



TF 



B. Low energies and effective spacetime metric 

The mode functions K of excitations with low en- 
ergies, ft n <C f-io, are localized inside the condensate. 
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Hence, changing the bounds of the integrals ([26[) and 
([271 from infinity to the Thomas-Fermi radius xtf will 
not alter these matrix elements significantly. Bearing fur- 
ther in mind that x 2 /2 + ,gof?o F — Mo = for x < xtf, it 
follows Afnm — and A4 nm — J d D x h+2goQoh^ n . Also 
A n 5 nm = J d D x h+2g g h^ n because of the restriction to 
low energies, fl n = V 'A n B n <C /io ~ goQo, and one gets 
M-nm = A n 5 nm . Hence, the evolution equations of dif- 
ferent eigenmodes approximately decouple and I obtain 
second-order equations of motion for phase 



dr 2 



2—^— + f 2 A n B n 

OT OT 



0. 



and density fluctuations 

d 2 

dr 2 J 



0, 



(28) 



(29) 



where I also introduced proper time dr = dt/b 2 . Equa- 
tion l|28p is the evolution equation of a mode of a 
minimally-coupled massless scalar field in a Friedman- 
Lemaitre- Robertson- Walker spacetime [U, [22] , provided 
the scale factor oflrw of the space-time is identified with 
l/f, cf. 



a FLRW — J 



(30) 



[Note that the prefactor of the damping term is 2 in Eq. 
([28[l . while it is D in a D + 1-dimensional Friedman- 
Lemaitre- Robertson- Walker spacetime.] 

The analogy ([28[) is not restricted to the evolution 
equations in mode expansion but applies in the low- 
energy limit of the field equations lfTT[l as well: within 
the Thomas-Fermi approximation, x 2 /2 + g g — /i = 
for x < xtf and Vq = 0, and for low excitations, 
V^x+ ^ 4<7o0oX+j the phase fluctuations obey a second- 
order field equation [13] 



8t 2 
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= 0. 



(31) 



which is similar to that of a minimally coupled scalar field 
in a Friedman-Lemaitre-Robertson- Walker spacetime. 

Having established this kinematical analogy, cf. (|3l)l . 
some of the concepts of general relativity can be applied 
to time-dependent Bose-Einstein condensates. Sonic 
analogs of horizons [U, [!]], SI SHI l45l . l46l | are of par- 
ticular interest for the study of non-equilibrium effects, 
because they give a rough estimate whether and when 
adiabaticity will be violated and the (quantum) fluctua- 
tions freeze and get amplified, i.e., (quasi-)particle pro- 
duction occurs. An effective particle horizon occurs, if 
a phonon emitted at a time tq can only travel a finite 
(co-moving) distance, i.e., if the integral 



A = / c s (r')dr' = / f(r')dr' (32) 



converges to a finite value Anorizon for r(i — > oo). 
Wavepackets emitted at time r at the origin x = can 
reach only points within the horizon, x < Anorizon, in a fi- 
nite time. All other points are concealed by the horizon. 
[For simplicity, I assumed in Eq. (|32p an homogeneous 
sound velocity c s .] 



C. Particle production in static traps 

In order to point out the analogy of phase fluctua- 
tions to cosmic quantum fields, I formulated the evolu- 
tion equations ([28]) and l[3"Tj) using proper time r. On 
the other hand, experiments are usually performed in 
the laboratory and thus the variations of trap frequency 
u> and coupling strength g are prescribed in laboratory 
time t. Since r is a complicated function of i, it is not 
quite obvious whether or not the quantum fluctuations 
will experience non-adiabatic evolution for a given mod- 
ulation of uj(t) or g{t). In laboratory time t, Eq. |29|) 
reads 



dt 2 



Ml 

J b(t) dt 



Sg n = 0, 



(33) 



which is the evolution equation of a damped harmonic 
oscillator with time-dependent coefficients 2b/b and 



(34) 



r 



To 



Initially, when b = 0, the field modes perform free os- 
cillations. Upon the gradual increase of the damping 
term 2b /b with respect to the oscillation frequencies fl n , 
the non-adiabatic evolution of the quantum fluctuations 
slowly sets in, until they finally freeze and get amplified 
when both terms, 2b/b and fl n are of the same order 

Let me discuss the two extremal ways a time- 
dependent scale factor b(t) can be achieved: firstly, only 
the trap frequency might be varied, while the interac- 
tion strength g = 1. With instantaneous frequencies 
fi„ oc b~ x ~ D l 2 , adiabaticity can be violated for any fi- 
nite change b/b ^ if b becomes sufficiently large. The 
quantum fluctuations cannot adapt to the changing back- 
ground any more, they freeze and get amplified. And, 
secondly, for static traps, w = 1, where g(t) is time- 
dependent. Then, the situations is not so clear because 
fll oc g/b 2+D = 1+6/6, cf. Eq. J3j). Hence, only the rapid 
acceleration of the scale factor |6/6| > 0(1) will lead to 
notable changes of the excitation frequencies. On the 
other hand, adiabaticity could be violated by increasing 
the magnitude of the damping term, 26/6. Then, how- 
ever, a continuous acceleration of 6 is required, because 
otherwise, if 6/6 was constant, the system would equilib- 
riate. 

As an example for the absence of particle production 
inside a static trap, uj = 1,1 will consider an exponential 



6 



sweep of the coupling coefficient 

g(t) = exp{ 7 f} 



(35) 



with 7 > 0. In this case, Eq. |(3j) for the scale factor 
permits an analytic solution 



b(t) 



(2 + D) 



7 2 + (2 + D) 2 



exp 



It 



2 + D 



(36) 



The coefficients of Eq. (|33|) become time-independent 



20 



D 



(2 + D) 2 



A n B n , 



(37) 



and the density eigenmodes are just damped harmonic 
oscillators with solutions 



(38) 



where <5^ is some residual oscillating function with fre- 
quency Q,,,. Hence, the density-density fluctuations di- 
minish ((X^. t ) 2 ) cx e~ 2lt ^ 2+D \ and, consequently, the 
phase-phase fluctuations increase. But this is just the 
adiabatic evolution, because A n - t cx e 2lt '^ 2+D ' and thus 
(X+ t ) ad = Q n . t /(2A n . t ) cx e -^t/(2+D)^ rf App rgj Thig 
means that no quasi-particle production occurs for the 
dynamics ([351 in the hydrodynamic regime, i.e., for low- 
energy excitations with order parameter treated in the 
Thomas-Fermi approximation. These findings can also 
be inferred from the (absence of an) effective particle 
horizon (|32)l . For the particular shape <[36|l of the scale 
factor b follows 



A cx / /(f) 



dt 



dt 



oo . 



(39) 



because /(f) oc exp{27f/(2 + D)} cx b 2 {t). 

Note that the presented solution assumes g = e 7 * at 
all times, especially also for t < fj n - Hence g(t') ^ 1 
at some time t' when b(t') = 1, cf. Eqs. ([35]) and ([36| . 
On the other hand, b(t- m ) — 1 and <?(fi n ) = 1 for a con- 
densate at rest, see Eq. ([3]). Since both solutions (static 
initial state and exponential sweep) cannot be matched 
at fjn such that g and b are both continuous, the switch- 
ing on of the exponential sweep would excite breathing 
oscillations. These oscillations, however, generally affect 
particle production, e.g., through parametric resonance. 



V. QUASI-ONE-DIMENSIONAL CONDENSATE 

The simplest application of the presented formalism 
consists in a quasi-one-dimensional condensate. In highly 
anisotropic traps, where the perpendicular trap fre- 
quency uj± is much larger than the chemical potential, 
the motion in the perpendicular directions is restricted to 



the ground state and might be integrated out. An effec- 
tively one-dimensional field equation ([T]) follows, where 
the interaction strength 



Sid 



,93D 

2-na 2 



= 2a, ui i 



(40) 



can be varied through Feshbach resonance or by changing 
oj±. However, one should bear in mind that the transver- 
sal extent of the condensate a±_ = 1 / ^/moj± has to be 
much larger than the s-wave scattering length a s such 
that the interaction of different atoms can still be de- 
scribed through three-dimensional scattering theory. For 
simplicity, I will adopt in this section the Thomas-Fermi 
approximation ([9]) for the order parameter but will per- 
mit arbitrary energies for the excitations. 



A. Spectrum 
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FIG. 1: Frequencies of odd (blue) and even (red) excitations 
versus the chemical potential /^o calculated using 100 har- 
monic oscillator basis functions. The frequencies equal the 
chemical potential at the intersection with the black line. 
The lowest odd mode converges for fio —* oo to the classical 
Thomas- Fermi breathing frequency fibi^ath — %/3wo (green), 
which follows from Eq. $3§. 



In Figure [H the excitation frequencies fl n = V A n B n 
of the lowest modes are plotted versus the chemical po- 
tential ^o- For fiQ = 0, one has the equidistant spectrum 
of the harmonic oscillator Q n = (n + l/2)w - Whereas 
for high fiQ ^> ft n , the frequencies become almost in- 
dependent of the chemical potential. In particular, the 
frequency Sli of the lowest excitation with odd parity 
tends for [iq — > oo towards the Thomas-Fermi breathing 
frequency ^breath = "V^^o obtained from Eq. ([3]). 

The discrepancy between these two frequencies fl\ and 
^breath ^ 0T finite values of the chemical potential /j-o hints 
at shortcomings of the Thomas-Fermi approximation. 
In particular, Eq. |(3|) does not describe the breathing 
motion of the background properly and care must be 
taken when employing Eqs. ([9]) for the order parame- 
ter. Nonetheless, Eq. still predicts the correct order 
of magnitude of the characteristic response time of the 



7 



background, 1/wbrcath = ^(V w breath)j to variations of 
trap frequency or interaction strength. 

Hence, it is still possible to discuss several cases, where 
the difference of Wbrcath and ^breath i s either small or does 
not matter: firstly, if the shape of the condensate varies 
only slowly, i.e., if b <C Wbrcath and the condensate can 
adapt to changes of ui and/or g immediately. Secondly, 
if no breathing oscillations are excited, e.g., because the 
condensate expands or contracts, see also Eqs. f35|) and 
(|36ll . And, thirdly, for very large chemical potentials, 
fio — » oo, the Thomas-Fermi approximation becomes ex- 
act. The quantum fluctuations are in the hydrodynamic 
regime and their effective evolution equations i|28p de- 
couple. In order to address cases where breathing of the 
background occurs, it would be necessary to abandon the 
Thomas-Fermi approximation ^ and to solve the Gross- 
Pitaevskii equation ([6]), which could, e.g., be done by ex- 
panding the order parameter -00 into oscillator functions 



B. Exponential sweep in stationary condensate 

A stationary condensate, i.e., 6=1, can be accom- 
plished through simultaneous variation of trap frequency 
ui and interaction strength g, cf. Eq. 



fit) 



£(*) 
50 



(41) 



though, one should be aware that this only holds within 
the Thomas-Fermi approximation: the instantaneous 
chemical potential must at all times be much larger than 
the trap frequency 



/WW = Mo/ 2 (i) = wo/(i) 



(42) 



If both were of the same order, the kinetic term, 
— in the Gross-Pitaevskii equation |(6]) becomes 
important; no stationary background could be realized 
even for simultaneous variation of uj and g. 



Analytical effective spacetime solution 



the time-dependence l(43j) . the integral 11321) is finite and 
an effective sonic horizon occurs; the quantum fluctua- 
tions freeze and get amplified. 

Instead of solving Eq. lj44j) for 6<j), I will consider the 
evolution equation of the density fluctuations X+ 



X+ = 0, 



(45) 



where A n B n — j 2 e 2 " ttrl . Obviously, all modes undergo 
the same evolution just at different times. Eq. (I45J1 can 
be solved analytically in terms of Bessel functions [48] 



x: 



^{a„^0 + 4<V)}- (46) 



where z = —j(t — t n ). The Hankel functions Hq' 2 ^ have 
the proper asymptotics for early times t — > — oo such that 
the operators a n annihilate the initial vacuum state. The 
phase fluctuations X~ = (l/2B n )dX+ / dt read 



X „ 



^^{a^VR^ffjV)} • (47) 



From these expressions f46|) and l(47jl . I can infer the 
correlations of each mode. At late times follows 



(X+)^(t-oo) 



2-yB„ 

7T 

7 
2ttB v 



(t - t n f 



(48) 



Comparison with the adiabatic values ((X^) 2 ) a( j, cf. Eq. 
(|B8[) . yields the quasi-particle number at late times 



N n (t ^oo) = I e T(*-*»)-i. 



(49) 



The occupation number of all modes grows exponentially 
though at different times t — t n , where the shift t n is de- 
termined by the excitation frequencies V 'A n B n = je ltn . 



However, one should bear in mind that Eq. l(45|) is only 
valid for a limited time before leaving the hydrodynamic 
regime. 



For an exponential sweep 



2. Numerical results 



/ 2 (t) = e- 2 ^, 7 >0, (43) 



the effective second-order equation (|3T|l for the phase 
fluctuations in the hydrodynamic limit 



dt 2 



d 



0. 



(44) 



is that of a massless scalar field in a de Sitter spacetime 
with exponentially growing scale factor Qflrw = 1// = 
e 7 *, cf. Eq. 1(30]) - which is believed to describe the uni- 
verse during the epoch of cosmic inflation (2lL l4l| . For 



In order to go beyond the effective space-time de- 
scription and thus the analytical findings (|46ll - lf4"9j) . I 
will now consider the full evolution equations <[25|l . The 
sweep rate 7 = 0.1 shall be chosen such that all modes 
evolve adiabatically at first, i.e., 7 <§; £l n {tm) for all n. 
When subsequently reducing trap potential and coupling 
strength, the excitation frequencies Q n (t) decrease and 
non-adiabatic evolution sets in at different times for each 
mode. The fluctuations freeze and get amplified. 

Figure [2] shows the instantaneous particle numbers 
of the lowest three modes for initial chemical potential 
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fio = 50 and 7 = 0.1. The lowest excitation, n = 0, which 
becomes non-adiabatic first, acquires the largest particle 
number. The next two modes, n = 1,2, experience less 
squeezing, though, remarkably, N2 > N\ - an unexpected 
result, which can be explained by the coupling of differ- 
ent modes: the second even mode, n = 2, gets populated 
from the principle excitation, n = 0, whereas the cou- 
pling matrix elements between n = and n — 1 are zero 
because of different parity. 
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FIG. 2: Instantaneous particle number in the lowest three 
quasi-particle modes for an exponential sweep l[43|l . The ini- 
tial chemical potential is /io = 50 and the sweep rate 7 = 0.1. 
At t » 80, the chemical potential equals the trap frequency. 
For the calculation of the eigenmodes h n and the coupling 
matrices (| 1 8 j> . I used the harmonic oscillator functions h a (x) 
up to a = 79, cf. Appendix lIIII The lowest 20 modes X n were 
then propagated. The numerical accuracy was set to 10 -6 . 



the evolution of the quantum fluctuations non-adiabatic, 
since the expansion/contraction of the background might 
compensate for changes of g such that the sound velocity 
(in comoving coordinates) remains constant. Breathing 
oscillations of the background, excited, e.g., by the sud- 
den change of the interaction coefficient g, might still 
yield a notable amount of quasi-particles. And secondly, 
if the excitation energy is of the same order as the chem- 
ical potential, different quasi-particle modes couple. 

The amplification and freezing of the fluctuations and 
also the coupling of different modes was illustrated in 
an example, where the trap frequency and interaction 
strength were exponentially ramped down such that the 
shape of the condensate remains constant. For the con- 
sidered parameters, a quasi-particle number of 0.5 was 
obtained in the lowest mode, though higher occupation 
numbers could be achieved by faster sweep rates or start- 
ing with a higher chemical potential. The inversion of the 
occupation number in the next two modes could be at- 
tributed to the inter-mode coupling: although the third 
excitation, n = 2, experiences a much briefer period of 
non-adiabatic evolution than the second mode, n = 1, 
only the former couples to the lowest mode, n = 0, and 
gets populated from it. 
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VI. SUMMARY 

The main objective of this Article was the investigation 
of quantum fluctuations in time-dependent harmonically- 
trapped Bose-Einstein condensates with repulsive inter- 
actions. To this end, the linear fluctuations were ex- 
panded into their initial eigenmodes and the field equa- 
tions were diagonalized. This diagonal form, however, 
persists only as long as the condensate is at rest; as soon 
as trap frequency or interaction strength are varied, the 
coupling of different modes sets in. (Only part of which 
can be accounted for by transforming to the instanta- 
neous eigenmodes, though the definition of instantaneous 
eigenmodes is a non-trivial task.) 

Two regimes were identified: firstly, for energies much 
smaller than the chemical potential, the coupling of dif- 
ferent modes is negligible and an effective space-time 
metric might be introduced for the phase fluctuations. 
This, however, necessitates a redefinition of the time co- 
ordinate such that the required change of trap frequency 
and/or interaction strength for a certain dynamics of this 
effective space-time is not obvious. It turned out that 
the sole variation of the interaction coefficient g(t) in 
a smooth monotonic way is hardly sufficient to render 



APPENDIX A: EIGENMODES 

The aim of this Appendix is the derivation of the ini- 
tial eigenfunctions of density and phase fluctuations. 
To this end, let me expand the field operator into any or- 
thonormal basis {h a (x)} of the underlying Hilbert space 
L 2 (R D ) with some operator- valued coefficients xt{t) 

X ±{x,t) = h a {x)xt(t) ■ (Al) 

Lower-case Greek indices (a, /?,...) shall denote the com- 
ponents in this arbitrary basis, while lower-case Latin 
indices (to, n,...) label the initial eigenmodes. 

If the condensate is initially at rest, the scale factor 
b = 1 and the phase of the order parameter is homoge- 
nous V x 0o = such that vo = and thus V nm = 0, cf. 
(fl9l) . The evolution equations l[TTj) simplify considerably 
and can be expanded into the basis {h a }. It follows, cf. 

m 

d f 

_2 ^Xq = AxpXp = j d D x hJC+hpxjii ) 

Id f 

2di*" = Baf3 *P = / dDxh »fc~ h (3Xp , (A2) 
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where JC± are defined in Eq. ([12]) . The matrices A a p 
and B a p are real and symmetric but do not commute 
due to [£+,£_] ^ 0, cf. HU). From Eqs. {£2}, I obtain 
second-order evolution equations for phase and density 
fluctuations 



9 2 

d 2 „, , 



(A3) 



which can be diagonalized by a transformation to the 
eigenvectors. 

Because the matrix A a pBp 1 ^ B a pAp 1 is not symmet- 
ric, phase and density fluctuations obey different eigen- 
value equations 



B a pAp 7 v™ = \ n v2 ) 



(A4) 



where n denotes the vectors and a merely counts the 
components in the particular basis representation l|Aip . 
For brevity, I will only discuss the eigenvectors of AB in 
the following, though the same applies for those of BA 
as well. The eigenvectors are generally not orthogonal 



(A5) 



but {v™} usually still forms a basis. (Note, however, that 
the eigenvectors of a non-symmetric matrix not always 
span the entire vector space. But if the v n were no basis, 
the evolution equation (|A3[) could not be diagonalized. 
This would mean that there existed some fluctuations, 
which have constant losses to some eigenmodes - rather 
unphysical in view of the stationary initial state consid- 
ered here. Therefore, I will not discuss this case any fur- 
ther.) With the v n = h a v n a forming a basis of L 2 (R D ), 
there must exist a dual basis {v d > n } in the space of linear 
functionals (the dual) on L 2 (M. D ), which obeys 



i D xv d ' n v r 



(A6) 



[Roughly speaking, the elements of L 2 (R D ) become col- 
umn vectors in the basis expansion (|A1|) . whereas row 
vectors correspond to the functionals on L 2 (R D ), i.e., the 
elements of the dual. Since the eigenvectors v n form a 
basis, the matrix with the v n as columns must be invert- 
ible. The rows of the (left) inverse matrix then comprise 
the elements of the dual basis {v d ' n } in the particular 
basis expansion l|Al|) .] 

After the multiplication of the first of Eqs. (|A4|) with 



,,d,m 



from the left and summation over a follows 

V a ' AapBp^V^ = V a ' V a X n = nm A m , 



(A7) 



which, because {v n } is a basis, implies that the v - n are 
the left eigenvectors of A a pBp 1 with the same eigenvalues 
A n 



,,d,n 



A a pB 



01 



i> d >" \ 



(A8) 



Transposition yields 



An< 



(A9) 



i.e., the eigenvalue equation of BA, cf. (j A4f) . Hence, 
v n ex v d,n and AB and BA (i.e., density and phase fluctu- 
ations) must have the same spectrum {A„} = {A„}. For 



simplicity, v n 



. ,d,n 



, which can be achieved by renor- 



malization of v n , in the following. 

Note that the spectrum of a real non-symmetric matrix 
might contain pairs of complex conjugate eigenvalues A„, 
A* , which can be seen when taking the complex conju- 
gate of (| A4|) . Complex eigenvalues 3A„ ^ are associ- 
ated with exponentially growing solutions, i.e., unstable 
modes. Since I am interested in the quantization of the 
stationary initial state, I will not discuss this case but in- 
stead assume A„ el Vn. (As can be easily verified, the 
components of the eigenvectors t;™ and u™ must be real- 
valued as well.) Nonetheless, complex eigenvalues might 
still occur in dynamical situations, e.g., during the sig- 
nature change event proposed in [9j or during (quantum) 
phase transitions j4§] (see also Ref. [13] for an illustrative 
example) . 

The initial evolution equations (j A3[) can be diagonal- 
ized by multiplication with the left eigenvectors 



2- 

dt 2 

dt 2 



VaXa = "V^AapB^X-y = ~ Kv^X-y > 

-vZB af) AfrX+ = - A„v"x+ ■ (A10) 



V n Y + 



which leads to the definition of the density and phase 
fluctuation eigenmodes X„ through 



X- = y n $- x + = v n r + 

n ^aAa ) n ^aAa 

where the spatial mode functions 



(All) 



(A12) 



follow from comparison of Eqs. (|Al[) with l[T5|) and the 
duality (|A6[) of the u™ and w™ implies Eq. (fTS)) for the 

K • 

In view of (|Allj) . the initial evolution equations (|A2j) 
can be transformed to the new basis, cf. Eq. fTTj) 



d - _ - Id" 

_2 <9t^" = ^ nXn ' 2dt Xn = BnX ™ (Ai3) 

where the transformed matrices 



An m = V™A a pvJ = A n S„ 
Bnm. = V^BapV™ = B n S n 



(A14) 



have become diagonal. Note that the transformation 
(|Alip is not orthogonal because the matrices comprising 
of the eigenvectors w™ and are not orthogonal. Hence, 
the commutators are not preserved, in particular 



AapBf}-, 7^ BapAp-f , 
AnmBml — Bnm Ami 



(A15) 
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where the latter can be inferred from Eq. 1|A7|1 

(KS nm ) T = V^Bapvfj V^AjSvT = BnkAkm . (A16) 

Since A nm and B nm also commute with their product 
AnkBkm — A„5 nm , which is diagonal, they must be diag- 
onal, too. 

APPENDIX B: OBSERVABLES 



The duality condition l|T6|) does not fix the norm of the 
eigenvectors but still permits the multiplication with an 
arbitrary factor A n 



A n h 



+ 

n ' 



A„ 



(Bl) 



This renormalization then leads to a stretching/shrinking 
of the operator- valued coefficients, cf. lfT5|) 



x: 



1 x- 



x: 



Kx n 



(B2) 



Since this is merely a basis transformation, the time- 
evolution of the quantum fluctuations X„ must be un- 
affected. To see this, recall the definitions lfT8|) and 
(JH) of the coupling matrices: they are the matrix el- 
ements of the operators JC± and of vo"V x + (V x Vo)/2 
with respect to the basis functions and thus acquire 
additional factors as well. As expected, all of these 
factors cancel such that the time-evolution remains un- 
changed. In particu lar, the eigenf requencies are invariant 
fin = VA n B n -> y/AlA n B n /Al = fi n . 



1. Correlation functions 

Furthermore the observables should not depend on the 
particular normalization of the basis functions. The rel- 
ative density-density correlations at time t read 



(6g(x)Sg(x')) = h+ t (x)h+. t (x') 
(g{x))(g(x>)) y/go{x)Qo(x') 



y/g (x)g (x') 



(B3) 



where I used the instantaneous eigenmode basis h^. t , see 
Eq. lf22|) . Obviously, the factors A„ and 1 / A„ contributed 
by h^-t and X n . t cancel each other and the spatial correla- 
tions are independent of the normalization A„. Similarly, 
the expression for the spatial phase-phase correlations 



K;t( X ) h m ;t( X ') 

Veo(x)Qo(x') 
KA x ) h ™A x ') 

A^o( x )Qo( x ') 



(x n . t x n 



n\t m:t ' 



(B4) 



yields the same result regardless of the employed basis. 

Hence, the factors A n can be chosen at will. There ex- 
ist, however, several convenient choices for A„: firstly, 
the density modes might be normalized to unity, 
/ d D x(h+) 2 — 1. This is advantageous if left and right 
eigenvectors are the same, i.e., if AB is symmetric. The 
drawback is that if AB is not symmetric and therefore 
/i+ 7^ h~, only one of the mode functions, can be 
normalized to unity, whereas the norm of the h~ fol- 
lows from Eq. (JT6J) . And, secondly, these factors A„ 
can be fixed by demanding A n = B n = Q n . In this 
case, the prefactors of (|20|) initially obey F™(t- m ) — 
8nm/V2 and G™(t- m ) = iS nm /y2 and both quadratures 
((X^) 2 )(ti n ) — 1/2. However, one should note that nei- 
ther of the eigenfunctions is generally normalized to 
unity, fd D x(h±) 2 ?l. 

Of course, the coefficients (X^tX^-t) and (X~ t X^. t ) 
do depend on the particular choice of the basis functions 
h„. t and thus also on the factors A„. For instance, the 
adiabatic density-density and phase-phase correlations of 
a particular mode read 



X n , t X n , t ) (t) = cx -L , 

'ad 2il n . t Ai 



ad 



(*) 



fir. 



2A n: , 



■ocAf 



(B5) 



The dependence on the factor A„ becomes important 
regarding the low excitations in the Thomas-Fermi ap- 
proximation, cf. Sec. IIVI the modes decouple and it is 
not necessary to introduce a new spatial basis at the 
time of measurement. One has instead A n -t = f 2 A n - to 
and B n -t = B n - to such that (X+ t X+ t ) a d oc / while 
(X~. t X~. t ) a< i cx 1/f, i.e., the phase and density fluc- 
tuations apparently increase or decrease even for adia- 
batic evolution. In view of this A„ ambiguity, the (ab- 
solute) density-density or phase-phase correlations pro- 
vide no adequate measure for the squeezing (i.e., non- 
adiabaticity) of a single mode. The Fourier transforms 
of Eqs. I|B3|) and (|B4[) on the other hand are independent 
of An but do not represent the excitation eigenmodes. 



Bogoliubov transformation and particle 
production 



As will be shown in the following, the (instantaneous) 
quasi-particle number measures the relative deviation of 
density and phase fluctuations from their adiabatic val- 
ues. In view of the different expansions lfl5|) and |22|) 
of density and phase fluctuations into their initial and 
adiabatic eigenfunctions, the corresponding creation and 
annihilation operators d n , a„ and b n . t , b n; t, respectively, 
can be transformed by virtue of a Bogoliubov transfor- 
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mation. For the annihilators b n -t follows in particular 



b n -.t. = 



A n :; 



2A n -t ' y 2il n t 



(B6) 



with the Bogoliubov coefficients a nm and I3 nm . The first 
line follows from inversion of Eq. (|24|) and the second 
line can be inferred after transforming (f2Q|l to the new 



basis h n . t . Since (3 nm ^ for non-adiabatic evolution, 

the quasi-particle number operator N n (t) = b n . t b n; t will 
have a non-zero expectation value as well 



N n (t) 



t>n;tbn;t 



2A n , t 
i 

+ 2 



£|/W*)I 
I a 

Afl-.t / , 2 



- \2 



TO) 



(B7) 



where the commutator 



Noting that 



the prefactors Q, n -t/2A n; t and A n]t /2n n . t are just the adi- 
abatic density-density and phase-phase correlations, see 
Eq. l|B5p . the particle number can be rewritten 



N n (t) 



(TO) 2 } 



2Q~ 



(TO) 2 } 



2A„ 



2A n;t /n 

n:t 



2fln;t/A 



(B8) 



i.e., the instantaneous particle number gives just the rela- 
tive deviation of the density and phase correlations from 
their adiabatic values. Note that expression (|B8|) does 
not contain the correlations between different modes. To 
this end, it would be necessary to evaluate (N n N m ), 
which is fourth order in the b n -t. This observable can 
be reduced to expectation values quadratic in in the b n -,t 
by virtue of Wick's theorem, see, e.g., [5l| . 



3. Scaling 



Another interesting aspect regards the scaling of the 
correlation functions (|B3[) and (|B4[) with the interac- 
tion strength: within the Thomas-Fermi approximation, 
see Sec. llVl the evolution equations ([25| are indepen- 
dent of go- All properties of the linear excitations are 
determined by the chemical potential /iq and the vari- 
ations git) /go and ui(t)/u>o, in particular the expecta- 
tion values (X n . t X m . t ). Only the normalization fac- 
tor y/g (x)Q (x') = (l/g )y/[fi - V(x)][fio ~ V(x')} in 
Eqs. (|B3(1 and (|B4[) depends on go. Hence, the rela- 
tive density-density and phase-phase correlations in the 
Thomas-Fermi approximation are both proportional to 
go for fixed fiQ, 
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